First-passage times for random walks in bounded domains 
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We present a novel computational method of first-passage times between a starting site and a 
target site of regular bounded lattices. We derive accurate expressions for all the moments of this 
first-passage time, validated by numerical simulations. Their range of validity is discussed. We 
also consider the case of a starting site and two targets. In addition, we present the extension to 
continuous Brownian motion. These results are of great relevance to any system involving diffusion 
in confined media. 
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How long does it take for a drunkard to go from a 
given bar to another one? This time is known in the 
random walk literature as a first passage time (FPT), 
and it has generated a considerable amount of work for 
many years P, Q . The importance of FPT relies on the 
fact that many physical properties, including fluorescence 
quenching 3] , neuron dynamics 4] or resonant activation 
5] to name a few, are controlled by first passage events. 
Unfortunately, explicit determinations of FPT are most 
of the time limited to very artificial geometries, such as 
ID and spherically symmetric problems 

The determination of FPT for random walks in real- 
istic geometries is not just a theoretical challenge in its 
own right. It is actually a very general issue involved as 
soon as molecules diffuse in confined media as, for ex- 
ample, biomolecules diffusing in the cell and undergoing 
a series of transformations at precise regions of the cell. 
An estimation of the time needed to go from one point 
to another is then an essential step in the understanding 
of the kinetics of the whole process. 

Very recently, two important advances in the calcula- 
tion of FPT have been performed. First, in the case of 
discrete random walks, an expression for the mean first 
passage time (MFPT) between two nodes of a general 
network has been found |H. So far, however, no quan- 
titative estimation of the MFPT has been derived from 
this formula. Second, the leading behavior of MFPT of 
a continuous Brownian motion at a small absorbing win- 
dow of a general reflecting bounded domain has been 
obtained 0, 13 • In the case when this window is a small 
sphere within the domain, the behavior of MFPT has re- 
cently been derived 9J. This result is rigorous, but does 
not give access to the dependence of the MFPT with the 
starting site. 

In this letter, we present a new computational method 
that allows us to quantitatively extend all these results 
in three directions: (i) we obtain an accurate explicit for- 
mula for the MFPT, (ii) we also examine all the moments 
of the FPT, and (iii) we consider the case with two tar- 
gets. The method is presented in detail in the case of 
discrete random walks on regular lattices, and then the 
extension to continuous Brownian motion is outlined. 

We first consider a random walker on a bounded lat- 
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FIG. I: Modifications of the original lattice: arrows denote 
one-way links. 



tice, and we address the question of determining the 
mean time needed to reach one point of the lattice (target 
site T) from another one (starting site S). The bound- 
aries are assumed to be reflecting. The starting point of 
the method is a result known in the mathematical liter- 
ature as Kac's formula 10] . Indeed, after our previous 
work about first return times (FRTs) for random walks 
[Tll |. we found out that Kac's formula allows one to ex- 
tend our results to general finite graphs. Kac's result 
concerns irreductible graphs (ie from any point one can 
go to any other point), which admit a stationary proba- 
bility 7r(r) to be at site r. Let us consider random walks 
starting from a random point of a subset S of sites of the 
lattice, with a probability proportional to the stationary 
probability. Then, the mean FRT of the random walk, 
i.e. the mean number of steps needed to return to any 
point of S, is, according to Kac's formula, l/7r(S]), where 

This formula gives FRTs and not FPTs. However, we 
can use it to derive the MFPT (T) by slightly modifying 
the original lattice (see fig^l: we suppress all the original 
links starting from the target site T, and add a new one- 
way link from T to the starting point S. In this modified 
lattice, the FRT to T is just the FPT from 5 to T in 
the former lattice, plus one. In what follows, for sake 
of simplicity, we only consider regular 2D or 3D lattices, 
although the argument may be easily extended to any 
kind of graph. Let tt be the position of the target site. 
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and rs be the position of the starting site. Denoting 
7r(rT) — j, Kac's formula gives (T) = 1/j — 1. Thus, 
all we need to know is thus the stationary probability for 
the modified graph. It satisfies the following equation: 

^(r) = i ^ 7r(r') + J<5rr, - 

{r',r) (r'.rT) 

where (r,r') means that these two sites are neighbors 
and a is the number of nearest neighbors of a site (by 
convention, the sites on the boundaries are their own 
neighbors). The last two terms of the rhs of iQJ are due 
to the modifications of the lattice. To solve this equation, 
we define the auxiliary function tt', equal to tt for r ^ ry, 
with 7r'(r7-) = 0. It satisfies: 

n'{r) = -Y, 7r'iv')+jSrrs~jSrrr, (2) 

SO that it' has the following expression: 

7T'{r) = ^+jH{r\rs)-jH{r\rT), (3) 

where iV is the total number of sites, and H is the discrete 
pseudo-Green function 12], which is symmetrical in its 
arguments and satisfies: 

H{r\r')^^ ^ H(r"|r') + <5rr'-^ (4) 

(r",r) 

Indeed the solution satisfies equation ^ , and ensures 
that TT is a probability function (of sum unity). The con- 
dition 7r'(rT) = allows us to compute j and to deduce 
the following exact expression: 

(T) = N[H{rT\rT) - H{rT\rs)] (5) 

This formula is equivalent to the one found in 6] , but is 
expressed in terms of pseudo-Green functions. One ad- 
vantage of our method is that it may be easily extended 
to more complex situations, as we will show. Another ad- 
vantage is that, although the pseudo-Green function H is 
not known in general, it is well suited to approximations. 
The simplest one is to approximate the pseudo-Green 
function by its infinite-space limit, the "usual" Green 
function: H{r\r') ~ G{r — r'), which satisfies: 

G(r) = -Y^ G(r') + Sor- (6) 

(r',r) 

The value of G(0) and the asymptotic behaviour of G 
are well-known [l^. For instance, for the 3D cubic lat- 
tice, we have: G{0) = 1.516386 and G{r) ~ 3/(27rr) for r 
large. For the 2D square lattice, we have G(0) — G(r) ~ 
2/7r ln(r) + 3/7r In 2 -|- 27/7r, where 7 is the Euler gamma 
constant. These estimations of G are used for all the 
practical applications in the following. This infinite-space 
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FIG. 2: 3D - Influence of the distance between the source 
and the target. Simulations (red crosses) vs. theory (plain 
line). The domain is a cube of side 31, the target being in the 
middle of it. The source takes all the positions in a cube of 
side 15 centered on the target. 

approximation may be improved by two kinds of correc- 
tions. First, the constant term 1/A^ in equation Q may 
be taken into account: 

i7(r|r') c^G(r-r') + ^(r-r')^ (7) 

In the 3D case, this "uniform correction" is always weak: 
its order of magnitude is at most N^^^^. However, in the 
2D case, it is negligible only if 7Vln(r5— r^) S> (rg— r^)^. 

A second correction that may be taken into account 
is the influence of nearby boundaries. For flat bound- 
aries, it can be computed explicitly. Denoting by s(r) 
the symmetric point of site r with respect to the closest 
flat boundary, H becomes: 

H{r\r') ~ G(r - r') + G(r - s(r')) (8) 

If the boundary is not flat, this expression only gives the 
order of magnitude of the expected correction. These two 
alternative corrections correspond to two different ways 
to treat the effect of boundaries: {T)) is a mean-field type 
correction, whereas © is a local correction. One should 
use either (fTj) or ||SJ) mainly according to the position of 
the target. A rule of thumb, used in the following, is that 
as soon as one of the two corrections is negligible, the 
other one leads to good results. Indeed, the correction 
is useful for a target far from any boundary, whereas 
the correction ||SJ| is more appropriate if the target is 
close to a flat boundary. As for the limitations of these 
approximations, they are not to be used in two cases: (i) 
if neither nor ijHJl are negligible; (ii) if the target is 
close to an irregular boundary. 

We have compared the theoretical predictions with nu- 
merical simulations. We first checked (fig. |3) the behav- 
ior of the MFPT when the source-target distance varies 
(the FPT is averaged over 100 ODD random walks). We 
also studied the influence of the distance between the 
target and a boundary (figOll, using the correction ((SJ. 
Finally, we checked that our approximation was also cor- 
rect for the 2D case (fig^. Since in this case the uniform 
correction ([7|) is not negligible, we took it into account. 
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FIG. 3; 3D - Influence of the distance between the target 
and a flat boundary. Simulations (red crosses) vs. theory 
(plain line). The domain is a cube of side 41 whose center 
is at (0,0,0); the source is at (0,0,x-15) and the target is at 
(0,0,x-20). 



FIG. 5: 3D - Higher-order moments: Theory (black curve) 
vs. simulation (red crosses). The blue dotted curve is the 
moments of the exponential distribution whose average time 
is the MFPT. The n-th moment is normalized by N"; the 
domain is a cube of side 51 centered on the target at (0,0,0) 
and the source at (2,2,1). 
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FIG. 4: 2D - Influence of the distance between the source and 
the target. Simulations (red crosses) vs. theory: blue dotted 
line, without the uniform correction; black plain line (in the 
middle of the set of points) , with this correction. The domain 
is a square of side 61, the target being in the middle of it, and 
the source takes all possible positions. 



In all the cases studied, the numerical simulations val- 
idate our approximations. Thus, our method provides an 
efficient way for estimating the MFPT, which up to now 
was only known formally and for a few specific cases. 

Furthermore, it is possible to compute the higher-order 
moments of the FPT, using an extension of Kac's for- 
mula, which gives a relation between the Laplace trans- 
form of the FRT to a subset S, averaged on S, and the 
Laplace transform of the FPT to this same subset, aver- 
aged on the complementary subset S. 

.(E) ((e-T)^ - e-) = (1 - .(S)) (e- - l) (e^^)^ 

Both averages are weighted by the stationary probabil- 
ity TT. For any starting point AI different from T, the 
behaviour of the random walk from M to T is exactly 
the same on the original and modified lattices until it 
reaches T. Thus, the FPT from M to T is the same on 
both lattices. This remark allows one, by applying the 
above formula to the modified lattice and using the cor- 
respondance between the FRT to T and the FPT from 
S, to get a relation between the n-th moment of the FPT 



and the lower-order moments of the FPT: 

We denote by the stationary distribution of the mod- 
ified graph whose starting point is r. The lowercase r 
refers to the starting point of the walk. This allows one, 
by recurrence, to get an estimation of the n-th order mo- 
ment, for large enough domains, but in the 3D case only 
(in fact, H{r\r') has to be negligible when r' is far from 
r). After some calculations which will be detailed in a 
further publication, it can be shown that 

(T") = n!7V" [{Ho - HirslvT)) {Ho - i/)""' + 0{N-i) 

(9)' 

where_i7o = H{rT\rT) and H = l/iV^^ i/(r|r'). Note 
that H is independent of r' due to the symmetry prop- 
erty of H, and that H scales as N~3, since G(r) 1/r. 
A good estimation of H, to be used for practical appli- 
cations, is its value for a spherical domain, computed in 
the continuous limit, H = (18/5)(3/(47r))2/3Af-i/3. The 
estimations @ are confirmed by numerical simulations 

(fig. El . 

It should be pointed out that the moments ^ are close 
but not equal (see fig. Isj to the moments of an expo- 
nential distribution of the FPT. However, if the parti- 
cle starts randomly inside the volume, the moments are 
the same as those of the exponential distribution, with a 
correction proportional to N~'^^^ . This property sheds a 
new light on the quasi-chemical approximation , which 
assumes that, if a particle starts randomly in the volume, 
it has a constant exit probability at each time step, which 
leads to an exponential distribution. 

We now turn to the situation where the lattice contains 
several targets, relevant in many chemical applications 
0. For sake of simplicity, the calculation is driven in 
the case of two targets, but may be easily extended to 
more absorbing points. We compute here the eventual 
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Position X of tlie second target. 



FIG. 6: 2D: Two-target simulations. Simulations (red crosses) 
vs. theory (plain line). One target is fixed at (-5,0); The 
source is fixed at (5,0); The other target is at (x,3). The 
domain is a square of side 201, the middle is the point (0,0) 



hitting probability to a specified target Pi, the so-called 
" splitting probability" , as well as the mean time until 
the particle hits either of the two targets (T) . We modify 
the graph in the same way as in fig^ for both absorbing 
points, denoted by ri and r2, the bonds relating them 
to their neighbours become one-way bonds, and a link 
is added from each target to the starting point rg. We 
denote 7r(ri) = ji, 7r(r2) = j2 and j = ji + j2- Again, 
the relation (T) — 1/j — 1 provides the mean absorption 
time, and the probabilities to hit ri or r2, are respectively 
ji/j and j2/j- We obtain a relation analogous to 



^'(r) 



1-j 

N 



+ jH(r\rs)-jiH{r\r,)-j2H{r\r2), (10) 



then, writing 7r'(ri) = 7r'(r2) =0 



^ + {.n+j2)His-jiHoi-j2Hi2 = 

^ + {jl + h)H2s - ]2Ho2 - ]lHi2 - 



(11) 



where H12 — -ff(ri|r2) and, for i = 1 or 2, His — 
H{Yi\Ts), Hoi = H{ri\ri). These equations yield exact 
expressions for the mean absorption time and the split- 
ting probabilities, respectively: 



/r-j-i\ [Hoi —His][Hq2 — / 



H24-[Hl2-H2s][Hl2-Hl, 



Pi 



P2 = 



-"2s+-"01— -"12 



-Ho2 — 2H\2 



Again, these expressions give excellent results when com- 
pared to simulations (fig. 

We finally address the case of a continuous Brownian 
motion. The target T is now a sphere of radius a centered 
on Yt, and the Brownian particle has a diffusion coeffi- 
cient D. It still starts from the point rg, at a distance 
R from the center of the target. The results are quite 
similar to those obtained in the discrete case, and the 
details of the computation will be published in a future 
paper. The estimated MFPT within the infinite-space 
approximation are 



(Tsd) = 



V 



--4);(T2d) = 



(12) 



4:TtD \a 

where V and A are the volume and area of the domains. 
If the target is approximately centered, the uniform cor- 
rection gives a contribution to (T) of —R^/{6D) in 3D 
and —R^ /{AD) in 2D. The correction due to a fiat re- 
flecting boundary is the following: 



(Tsd) - 4^ (i + ^ - - if) 

(T2d) = 2^(lnf + lni^ 



(13) 



with d the distance between the center of the sphere 
T and the boundary, and R' the distance between this 
same center and the reflexion of the starting point by 
the boundary. Note that these results significantly ex- 
tend the (exact) formula of Pinsky |^, which only gives 
the leading term in a. 

In summary, we have presented here a new method 
of computation that yields very accurate expressions 
of mean first-passage times for discrete random walks 
and continuous Brownian motion. These approximations 
have proven to be especially useful when the target is 
roughly in the middle of the bounded domain or near 
to a fiat boundary. This approach also gives access to 
more complex quantities such as higher-order moments. 
These results may be of the greatest interest for systems 
involving diffusion in confined media. 
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